Información

El código y los archivos orginales de esta PEC, así como los informes en formato HTML pueden consultarse en el siguiente repositorio de GitHub. El repositorio es privado, así que se tendrá que requerir el acceso para poder ver el código.

Ejercicios

Ejercicio 1

En un estudio sobre una enfermedad del hígado se recogieron datos de 583 pacientes del departamento de digestivo de un hospital. De ellos, 416 experimentan una afección del hígado y 167 tienen otros problemas pero no se ven afectados por dicha afección. El archivo dataset.csv contiene un conjunto de variables que se consideran “marcadores” de dicha enfermedad, además del género de los pacientes. Las variables registradas son:

  • Age: Edad del paciente
  • Gender: Género del paciente
  • TB: Bilirubina total
  • DB: Bilirubina directa
  • Alkphos: Fosfatasa alcalina
  • Sgpt: Alamina Aminotransferasa
  • Sgot: Aspartato Aminotransferasa
  • TP: Proteínas totales
  • ALB: Albúmina
  • A/G: Ratio entre Albúmina y Globulina
  • Selector: campo utilizado para romper los datos en dos grupos etiquetados 1 (liver) o 2 (control) por un grupo de expertos.
data <- read.csv("dataset.csv") %>% dplyr::mutate(Selector = factor(Selector))
head(data)

1.1. Exploración de los datos

A)

Realizar un resumen numérico y gráfico de los datos. Observar si hay valores faltantes y tenerlo en cuenta en el análisis.

La manera más fácil de realizar un resumen numérico es usar la función summary().

summary(data)
##       Age           Gender          TB               DB        
##  Min.   : 4.00   Female:142   Min.   : 0.400   Min.   : 0.100  
##  1st Qu.:33.00   Male  :441   1st Qu.: 0.800   1st Qu.: 0.200  
##  Median :45.00                Median : 1.000   Median : 0.300  
##  Mean   :44.75                Mean   : 3.299   Mean   : 1.486  
##  3rd Qu.:58.00                3rd Qu.: 2.600   3rd Qu.: 1.300  
##  Max.   :90.00                Max.   :75.000   Max.   :19.700  
##                                                                
##     Alkphos            Sgpt              Sgot              TP       
##  Min.   :  63.0   Min.   :  10.00   Min.   :  10.0   Min.   :2.700  
##  1st Qu.: 175.5   1st Qu.:  23.00   1st Qu.:  25.0   1st Qu.:5.800  
##  Median : 208.0   Median :  35.00   Median :  42.0   Median :6.600  
##  Mean   : 290.6   Mean   :  80.71   Mean   : 109.9   Mean   :6.483  
##  3rd Qu.: 298.0   3rd Qu.:  60.50   3rd Qu.:  87.0   3rd Qu.:7.200  
##  Max.   :2110.0   Max.   :2000.00   Max.   :4929.0   Max.   :9.600  
##                                                                     
##       ALB             A.G         Selector
##  Min.   :0.900   Min.   :0.3000   1:416   
##  1st Qu.:2.600   1st Qu.:0.7000   2:167   
##  Median :3.100   Median :0.9300           
##  Mean   :3.142   Mean   :0.9471           
##  3rd Qu.:3.800   3rd Qu.:1.1000           
##  Max.   :5.500   Max.   :2.8000           
##                  NA's   :4

El resumen numérico nos puede dar varias informaciones, donde tenemos que la edad de los pacientes oscila entre 4 y 90 años, siendo la media y mediana muy similares (44.75 y 45), entre otras variables.

Buscando si hay valors faltantes (indicados como NA), se puede observar fácilmente que en el ratio de albúmina y globulina (variable A.G), hay un total de 4 valores que no están presentes.

Para proceder con el análisis, dado que el número de valores NA és muy peqeño y solo afecta a una variable, podríamos usar el método de complete-case analysis usando solo los casos que tengan todas las variables completas o el método de available-case analysis. Debido a su mayor senzillez, vamos a proceder con el filtraje de los casos incompletos, quedandonos con esos que tengan información para todas las variables.

data_filt <- data[which(complete.cases(data)),]
summary(data_filt)

Para hacer el resumen gráfico de las variables, podemos usar la función plot(), que nos dibujará una serie gráficos de dispersión relacionando las variables una contra una.

plot(data, main = "Scatterplot comparing all the variables one-vs-one")
Figura 1.1. Gráfico de dispersión que compara todas las variables (por pares) de nuestro set de datos. Los datos usados para este gráfico incluyen los casos incompletos y todas las variables, sean numéricas o no.

Figura 1.1. Gráfico de dispersión que compara todas las variables (por pares) de nuestro set de datos. Los datos usados para este gráfico incluyen los casos incompletos y todas las variables, sean numéricas o no.

A simple vista, las variables indicativas de la bilirubina directa y la bilirubina total (DB y TB, respectivamente) parecen tener una correlación bastante buena, así como las variables TP y ALB o ALB y A/G.




B)

Estudiar la relación entre las variables numéricas, especialmente entre la bilirubina total y la directa. ¿Están correlacionadas?

Para hacerlo sencillo, vamos a crear otro set de datos a partir de los datos filtrados (casos completos) con solo las variables numéricas (todas menos Gender y Selector)

data_num <- data_filt %>% dplyr::select(everything(), -Gender, -Selector)

En el resumen gráfico del apartado anterior, a simple vista parece que la bilirubina total (TB) y la directa (DB) tienen una buena correlación.

Seguidamente vamos a comprobar la correlación entre variables usando la función cor() y a dibujar un gráfico de correlación con la función corrplot().

corrplot::corrplot(corr = cor(data_num), type = "lower", 
                   number.digits = 2, method = "color", 
                   addgrid.col = "gray50", addCoef.col = "black",
                   title = "Correlation plot showing the correlation between each pair of variables")
Figura 1.2. Correlaciones entre las distintas variables. A más intensidad de color, más correlación, ya sea positiva (azul) o negativa (rojo). Para hacer este gráfico solo se han usado las variables numéricas y los casos completos.

Figura 1.2. Correlaciones entre las distintas variables. A más intensidad de color, más correlación, ya sea positiva (azul) o negativa (rojo). Para hacer este gráfico solo se han usado las variables numéricas y los casos completos.

En el gráfico anterior, donde a más intensidad de color más correlación, vemos que, en general, la correlación entre las distintas variables es pequeña (coeficiente de correlación de Pearson de menos de 0.3 en valor absoluto). Sin embargo, hay cuatro pares de variables que sí que tienen una alta correlación, siendo la bilirubina total i la bilirubina directa las variables más correlacionadas (Pearson = 0.87):

  • DB vs TB: Pearson = 0.87
  • Sgot vs Sgpt: Pearson = 0.79
  • ABL vs TP: Pearson = 0.78
  • A.G vs ALB: Pearson = 0.69

Mirando estos números (figura 1.2), además del gráfico de dispersión (figura 1.1), podemos decir que la bilirubina total y la directa están altamente correlacionadas.




C)

Hacer un gráfico multivariante de las correlaciones dos a dos.

Mirar apartado anterior




D)

Obtener una lista de las correlaciones dos a dos, ordenadas de mayor a menor.

En la Tabla 1.1 se muestran las correlaciones dos a dos, ordenadas de mayor a menor valor de Pearson en valor absoluto. Por razones de espacio, la Tabla 1.1 se encuentra cortada a 10 entradas y la lista completa puede verse en la Tabla S1.1 del material suplementario.

cor_list_1d <- cor(data_num) %>% 
  reshape2::melt() %>%
  dplyr::mutate(abs_pearson = sqrt(value^2)) %>%
  dplyr::arrange(desc(abs_pearson)) %>%
  magrittr::set_colnames(c("Variable 1", "Variable 2", "Pearson's correlation", "Absolute Pearson's"))

cor_list_1d %>% head(10) %>% knitr::kable(align = "c", caption = "Tabla 1.1. Lista de correlaciones dos a dos, ordenadas de mayor a menor correlación (en valor absoluto). Los valores de 1, indicando que una variable se enfrenta a ella misma, son mantenidos. Solo se muestran las 10 primeras entradas. La lista con todos los valores se puede encontrar en el material suplementario (Tabla S1.1).")
Tabla 1.1. Lista de correlaciones dos a dos, ordenadas de mayor a menor correlación (en valor absoluto). Los valores de 1, indicando que una variable se enfrenta a ella misma, son mantenidos. Solo se muestran las 10 primeras entradas. La lista con todos los valores se puede encontrar en el material suplementario (Tabla S1.1).
Variable 1 Variable 2 Pearson’s correlation Absolute Pearson’s
Age Age 1.000000 1.000000
TB TB 1.000000 1.000000
DB DB 1.000000 1.000000
Alkphos Alkphos 1.000000 1.000000
Sgpt Sgpt 1.000000 1.000000
Sgot Sgot 1.000000 1.000000
TP TP 1.000000 1.000000
ALB ALB 1.000000 1.000000
A.G A.G 1.000000 1.000000
DB TB 0.874481 0.874481




E)

Repetir el punto anterior considerando las correlaciones parciales entre pares de variables cuando se fijan todas las demás. ¿Se ven afectadas vuestras conclusiones?

Para calcular las correlaciones parciales, usaremos la función pcor() del paquete ppcor().

Primeramente haremos un gráfico de correlación como el de la Figura 1.2, con tal de observar las diferencias. Este gráfico podrá ser observado en la Figura S1.1 en el material suplementario.

par_cor_1e <- ppcor::pcor(data_num)

Viendo la Figura S1.1, se ve claramente que las correlaciones parciales cambian comparadas con las correlaciones obtenidas en la Fig. 1.2. Por ejemplo vemos como la correlación entre la bilirubina total (TB) y la directa (DB) disminuye ligeramente, mientras que las correlaciones parciales ALB vs TP (Pearson’s = 0.82), ALB vs A.G (Pearson’s = 0.90), y TP vs A.G (Pearson’s = -0.69), augmentan (en valor absoluto) mucho comparadas con las correlaciones anteriores.

Igual que en el apartado d, la lista completa de las correlaciones parciales se puede observar en la Tabla S1.2 en el material suplementario mientras que la Tabla 1.2 solo muestra las diez primeras entradas.

par_cor_1e_list <- par_cor_1e$estimate %>% 
  reshape2::melt() %>%
  dplyr::mutate(abs_pearson = sqrt(value^2)) %>%
  dplyr::arrange(desc(abs_pearson)) %>%
  magrittr::set_colnames(c("Variable 1", "Variable 2", "Pearson's correlation", "Absolute Pearson's"))

par_cor_1e_list %>% head(10) %>% knitr::kable(align = "c", caption = "Tabla 1.2. Lista de correlaciones parciales dos a dos, ordenadas de mayor a menor correlación (en valor absoluto). Los valores de 1, indicando que una variable se enfrenta a ella misma, son mantenidos. Solo se muestran las 10 primeras entradas. La lista con todos los valores se puede encontrar en el material suplementario (Tabla S1.2).")
Tabla 1.2. Lista de correlaciones parciales dos a dos, ordenadas de mayor a menor correlación (en valor absoluto). Los valores de 1, indicando que una variable se enfrenta a ella misma, son mantenidos. Solo se muestran las 10 primeras entradas. La lista con todos los valores se puede encontrar en el material suplementario (Tabla S1.2).
Variable 1 Variable 2 Pearson’s correlation Absolute Pearson’s
Age Age 1.0000000 1.0000000
TB TB 1.0000000 1.0000000
DB DB 1.0000000 1.0000000
Alkphos Alkphos 1.0000000 1.0000000
Sgpt Sgpt 1.0000000 1.0000000
Sgot Sgot 1.0000000 1.0000000
TP TP 1.0000000 1.0000000
ALB ALB 1.0000000 1.0000000
A.G A.G 1.0000000 1.0000000
ALB TP 0.8971925 0.8971925




F)

Calcular el vector de medias, la matriz de covarianzas y la matriz de correlaciones de las variables respuesta numéricas conjuntamente pero por separado para los dos grupos (niveles del factor Selector sin tener en cuenta el género).

Primero separamos los datos según el valor en la variable Selector (1 o 2).

data_num_selector1 <- data_filt %>% dplyr::filter(Selector == 1) %>% dplyr::select(everything(), -Gender, -Selector)
data_num_selector2 <- data_filt %>% dplyr::filter(Selector == 2) %>% dplyr::select(everything(), -Gender, -Selector)


mean_vect_selector1 <- colMeans(data_num_selector1)
mean_vect_selector2 <- colMeans(data_num_selector2)

cov_mat_selector1 <- cov(data_num_selector1)
cov_mat_selector2 <- cov(data_num_selector2)

corr_mat_selector1 <- cor(data_num_selector1)
corr_mat_selector2 <- cor(data_num_selector2)

Después de haber separado los datos, calculamos los vectores de medias con colMeans(), la matriz de covarianzas con cov() o la matriz de correlaciones con cor():

Vector de medias:

  • Selector == 1 (“Liver”):

46.1449275, 4.1804348, 1.9316425, 319.5362319, 99.97343, 138.173913, 6.4586957, 3.0584541, 0.9141787

  • Selector == 2 (“Control”):

41.3636364, 1.1448485, 0.3963636, 220.6848485, 33.8363636, 40.7636364, 6.5393939, 3.3393939, 1.0295758

Matriz de covarianzas:

  • Selector == 1 (“Liver”):

246.1871776, -2.6722181, -1.7319577, 376.866407, -436.8823385, -243.4586799, -3.1780187, -3.4973541, -1.3308976, -2.6722181, 51.2423524, 19.9290225, 318.8957469, 279.9389462, 509.423992, 0.0609333, -1.2364565, -0.4557098, -1.7319577, 19.9290225, 10.3203765, 164.7544198, 137.6199711, 247.9961785, 0.0670728, -0.5629437, -0.1916362, 376.866407, 318.8957469, 164.7544198, 7.227736310^{4}, 5390.7891006, 1.315774210^{4}, -5.9794926, -29.9190722, -18.3494859, -436.8823385, 279.9389462, 137.6199711, 5390.7891006, 4.546146410^{4}, 5.675900710^{4}, -10.6640541, -0.8768935, 1.9592638, -243.4586799, 509.423992, 247.9961785, 1.315774210^{4}, 5.675900710^{4}, 1.143361110^{5}, -7.4678598, -18.6288346, -5.948743, -3.1780187, 0.0609333, 0.0670728, -5.9794926, -10.6640541, -7.4678598, 1.2040283, 0.6576987, 0.0720592, -3.4973541, -1.2364565, -0.5629437, -29.9190722, -0.8768935, -18.6288346, 0.6576987, 0.6200131, 0.170283, -1.3308976, -0.4557098, -0.1916362, -18.3494859, 1.9592638, -5.948743, 0.0720592, 0.170283, 0.1063755

  • Selector == 2 (“Control”):

291.0133038, 0.2049335, 0.1732816, -190.6956763, -46.7998891, -61.6452328, -3.2686807, -2.2266075, -0.2056375, 0.2049335, 1.0193178, 0.5197982, 80.9239763, 8.5153104, 14.2240798, -0.1654361, -0.1452531, -0.0472918, 0.1732816, 0.5197982, 0.2724257, 41.6836031, 4.3908647, 7.3869401, -0.079429, -0.0732705, -0.0248308, -190.6956763, 80.9239763, 41.6836031, 2.00301210^{4}, 1341.8261641, 1384.5409091, -4.3997044, -16.1198263, -9.8257443, -46.7998891, 8.5153104, 4.3908647, 1341.8261641, 632.332816, 610.3452328, 0.9814856, 0.8766075, 0.0663326, -61.6452328, 14.2240798, 7.3869401, 1384.5409091, 610.3452328, 1336.8645233, -4.0711197, -2.3125831, 0.2007528, -3.2686807, -0.1654361, -0.079429, -4.3997044, 0.9814856, -4.0711197, 1.1094752, 0.7056338, 0.0987973, -2.2266075, -0.1452531, -0.0732705, -16.1198263, 0.8766075, -2.3125831, 0.7056338, 0.6061826, 0.1649558, -0.2056375, -0.0472918, -0.0248308, -9.8257443, 0.0663326, 0.2007528, 0.0987973, 0.1649558, 0.0825138

Matriz de correlaciones:

  • Selector == 1 (“Liver”):

1, -0.0237917, -0.0343603, 0.0893416, -0.13059, -0.0458882, -0.1845888, -0.2830782, -0.2600706, -0.0237917, 1, 0.86661, 0.165704, 0.1834117, 0.2104617, 0.0077575, -0.2193632, -0.195188, -0.0343603, 0.86661, 1, 0.1907604, 0.2009148, 0.2282998, 0.0190275, -0.2225444, -0.182898, 0.0893416, 0.165704, 0.1907604, 1, 0.0940437, 0.14474, -0.0202696, -0.141334, -0.2092676, -0.13059, 0.1834117, 0.2009148, 0.0940437, 1, 0.7872658, -0.0455808, -0.0052231, 0.0281741, -0.0458882, 0.2104617, 0.2282998, 0.14474, 0.7872658, 1, -0.0201273, -0.069967, -0.0539402, -0.1845888, 0.0077575, 0.0190275, -0.0202696, -0.0455808, -0.0201273, 1, 0.7612165, 0.2013494, -0.2830782, -0.2193632, -0.2225444, -0.141334, -0.0052231, -0.069967, 0.7612165, 1, 0.6630558, -0.2600706, -0.195188, -0.182898, -0.2092676, 0.0281741, -0.0539402, 0.2013494, 0.6630558, 1

  • Selector == 2 (“Control”):

1, 0.0118988, 0.0194613, -0.0789846, -0.1090977, -0.0988324, -0.1819103, -0.167643, -0.0419645, 0.0118988, 1, 0.9864065, 0.5663444, 0.3354075, 0.3853237, -0.1555667, -0.184786, -0.1630677, 0.0194613, 0.9864065, 1, 0.5642862, 0.3345439, 0.3870765, -0.1444763, -0.1803032, -0.1656165, -0.0789846, 0.5663444, 0.5642862, 1, 0.377035, 0.2675595, -0.0295136, -0.1462908, -0.2416909, -0.1090977, 0.3354075, 0.3345439, 0.377035, 1, 0.6638332, 0.0370555, 0.0447745, 0.0091831, -0.0988324, 0.3853237, 0.3870765, 0.2675595, 0.6638332, 1, -0.1057089, -0.0812366, 0.0191141, -0.1819103, -0.1555667, -0.1444763, -0.0295136, 0.0370555, -0.1057089, 1, 0.8604365, 0.3265298, -0.167643, -0.184786, -0.1803032, -0.1462908, 0.0447745, -0.0812366, 0.8604365, 1, 0.7375689, -0.0419645, -0.1630677, -0.1656165, -0.2416909, 0.0091831, 0.0191141, 0.3265298, 0.7375689, 1




G)

Calcular la varianza generalizada, la varianza total y el coeficiente de dependencia global \(\eta^2 = 1−|R|\) para cada grupo por separado.

Varianza generalizada: se defene como el determinante de la matriz de varianzas. Así pues, con la función det() y la función var() vamos a calcular la varianza generalizada (ex. det(var(data))).

var_g_total     <- det(var(data_num))
var_g_selector1 <- det(var(data_num_selector1))
var_g_selector2 <- det(var(data_num_selector2))

La varianza generalizada del total de los datos es 5.98808510^{15}, mientras que para los datos etiquetados como Selector == 1 es 3.202195510^{16} y para los datos etiquetados con Selector == 2 es 1.053902710^{7}.

Varianza total: se define como la suma de las varianzas de cada variable o como la traza de la matriz de varianzas y covarianzas. Como las varianzas son la diagonal de la matriz de varianzas-covarianzas, con la suma de la diagonal vamos a poder calcular la varianza total.

vt_global    <- sum(diag(var(data_num)))
vt_selector1 <- sum(diag(var(data_num_selector1)))
vt_selector2 <- sum(diag(var(data_num_selector2)))

La varianza total del total de los datos es 1.772031510^{5}, mientras que para los datos etiquetados como Selector == 1 es 2.323846210^{5} y para los datos etiquetados con Selector == 2 es 2.22934210^{4}.

Coeficiente de dependencia global o \(\eta^2 = 1−|R|\): para calcular el coeficiente de dependencia global tenemos que calcular el determinante \(||\) de la matriz de correlaciones \(R\) (\(|R|\)) y restarlo de 1.

coef_dep           <- 1-det(cor(data_num))
coef_dep_selector1 <- 1-det(cor(data_num_selector1))
coef_dep_selector2 <- 1-det(cor(data_num_selector2))

El coeficiente de dependencia del total de los datos es 0.9941683, mientras que para los datos etiquetados como Selector == 1 es 0.9917558 y para los datos etiquetados con Selector == 2 es 0.9998612.




1.2. Visualización de los datos en dimensión reducida

A)

Realizar un análisis de componentes principales y estudiar la calidad de la representación por diversos criterios según el número de dimensiones.

Primeramente, debemos tener en cuenta que el análisis de componentes principales (en adelante PCA) puede hacerse a partir de la matriz de covarianzas o la matriz de correlaciones.

Normalmente, si las distintas variables están medidas en escalas muy diferente y la varianza de las variables difiere mucho una de otra, no es recomendable que se ejecute el PCA a partir de la matriz de covarianzas, pues las primeras componentes principales van a estar definidas en base a una o pocas variables con una gran varianza.

A continuación observamos las varianzas de cada variable (presentes en la diagonal de la matriz de covarianzas):

diag(var(data_num)) %>% round(3) %>% knitr::kable(caption = "Tabla 1.3. Varianzas de las distintas variables de nuestro set de datos.")
Tabla 1.3. Varianzas de las distintas variables de nuestro set de datos.
x
Age 263.146
TB 38.784
DB 7.933
Alkphos 59322.381
Sgpt 33555.955
Sgot 84013.042
TP 1.176
ALB 0.631
A.G 0.102

Como se puede observar en la Tabla 1.3, las varianzas de las distintas variables son muy diferentes, siendo para algunas variables (i.e. Alkphos) mucho mayor que para otras (i.e. TP), por lo tanto, lo más adecuado en este caso sería usar la matriz de correlaciones para hacer el PCA.

Para buscar las componentes principales, usaremos la función princomp(), a la qual debemos dar nuestro set de datos (con únicamente las variables numéricas sin NA) e indicar que queremos que haga el cálculo en base a la matriz de correlaciones (cor=T).

pca <- princomp(x = data_num, cor = T)

pca_vars <- pca$sdev**2

pca_sum <- summary(pca)

pca_sum
## Importance of components:
##                           Comp.1    Comp.2    Comp.3    Comp.4     Comp.5
## Standard deviation     1.6594914 1.4236792 1.1685377 0.9787747 0.91905411
## Proportion of Variance 0.3059902 0.2252069 0.1517200 0.1064444 0.09385116
## Cumulative Proportion  0.3059902 0.5311971 0.6829171 0.7893616 0.88321272
##                            Comp.6     Comp.7     Comp.8      Comp.9
## Standard deviation     0.81632677 0.45105044 0.35470439 0.235445205
## Proportion of Variance 0.07404327 0.02260517 0.01397947 0.006159383
## Cumulative Proportion  0.95725598 0.97986115 0.99384062 1.000000000
prop_var_pca <- data.frame(pc = c("PC1", "PC2", "PC3", "PC4", "PC5", "PC6", "PC7", "PC8", "PC9"),
                           prop_var = c(0.3059902, 0.2252069, 0.1517200, 0.1064444, 0.09385116, 
                                        0.07404327, 0.02260517, 0.01397947, 0.006159383),
                           acum_var = c(0.3059902, 0.5311971, 0.6829171, 0.7893616, 0.88321272, 
                                        0.95725598, 0.97986115, 0.99384062, 1.000000000))
  
ggplot(prop_var_pca) + 
  
  geom_col(aes(x = pc, y = prop_var, group = 1, fill = pc), position = "dodge", show.legend = F) +
  geom_point(aes(x = pc, y = acum_var/3,), show.legend = F) +
  geom_line(aes(x = pc, y = acum_var/3, group = 1), show.legend = F) +
  
  scale_y_continuous(name = expression("Explained variance (proportion)"), 
                     sec.axis = sec_axis(~ . * 3 , name = "Cumulative variance (proportion)"), limits = c(0, 0.35)) + 
  
  xlab("Principal component") +
  ggtitle("Variance explained by each principal component")
Figura 1.3. Proporción de la varianza explicada por cada componente principal (barras, eje izquierdo) junto con la proporción acumulada (puntos y líneas, eje derecho). Este tipo de gráfico a veces se llama 'scree plot'

Figura 1.3. Proporción de la varianza explicada por cada componente principal (barras, eje izquierdo) junto con la proporción acumulada (puntos y líneas, eje derecho). Este tipo de gráfico a veces se llama ‘scree plot’

¿Cuántas componentes principales debemos escoger? Para escoger con cuantas componentes principales nos quedamos para continuar el análisis, podemos usar varios criterios:

  • Criterio de porcentaje: este criterio se basa en definir un porcentaje mínimo, por ejemplo 80%, de la varianza que queremos explicar con nuestras componentes principales:

    • En el caso de definir un porcentaje mínimo del 80%, nos quedaremos con las primeras 5 componentes principales, ya que así pasamos del 80% de varianza explicada, concretamente tenemos un 88.3% (Figura 1.3).

    • Sin embargo, en algunos casos, el punto de corte basado en el porcentaje se define cuando, al considerar una componente principal más, no se aumenta mucho en el porcentaje de varianza explicada. En otras palabras, hay una estabilización en el porcentaje de varianza explicada. Así pues, fijándonos en la Figura 1.3, se puede ver como de la componente principal 5 a la componente principal 6 aún hay un aumento considerable de la varianza explicada, mientras que a partir de la 6, la proporción de varianza acumulada aumenta muy poco con cada componente que añadimos, por lo que tendríamos razones para coger hasta 6 componentes principales.

  • Criterio de Kaiser: al usar la matriz de correlaciones para obtener las componentes principales, estamos asumiendo que las variables observables tienen varianza 1, por lo que una componente principal con una varianza por encima de 1 explica más variación de los datos que cualquier variable observable, mientras que una componente principal con varianza menor que 1 explica una menos variación que una de las variables observables y, por lo tanto, la descartaríamos para el análisis.

    • Siguiendo este criterio, escogeríamos las 3 primeras componentes principales (desviación estándard de 1.6594914, 1.4236792 y 1.1685377 respectivamente), ya que son las únicas con una desviación estándard mayor que 1 y, por ende, varianza mayor que 1. En este caso, con 3 componentes principales, explicaríamos un 68% de la varianza de los datos.

    • En algunos casos, se considera que el punto de corte debe ser 0.7, por lo que en este caso, nos quedaríamos con 5 componentes principales, dado que la 6a componente principal tiene una varianza de 0.6663894.

Calidad de la representación:

Para ver la calidad de la respresentación podemos usar la función get_pca_var() del paquete factoextra, que nos devuelve los resultados (coordenadas, cosenos cuadrados y contribuciones) de las variables en el PCA. En este caso nos interesan los cosenos cuadrados (cos2), que representan la calidad de la representación de las variables y es calculado como las coordenadas (coords) al cuadrado. En este caso, un valor alto de cos2 nos indica una buena representación de la variable en la dimensión o componente principal en cuestión.

pca_var <- factoextra::get_pca_var(res.pca = pca)

corrplot::corrplot(pca_var$cos2, is.corr=FALSE, method = "color", addgrid.col = "gray50", addCoef.col = "black",
                   title = "Cos2 of variables in each dimension or PC.")
Figura 1.4. Cos2 (calidad de la representación) de las variables en cada dimensión o componente principal. Un valor alto de 'cos2' indica que la variable tiene una buena representación en la dimensión en cuestión.

Figura 1.4. Cos2 (calidad de la representación) de las variables en cada dimensión o componente principal. Un valor alto de ‘cos2’ indica que la variable tiene una buena representación en la dimensión en cuestión.

En la Figura 1.4 podemos ver la calidad de la representación de cada variable en cada dimensión: a un valor más alto de cos2 mejor representación tiene dicha variable. Por ejemplo, en la primera dimensión (o componente principal), las variables mejor representadas son la TB, la DB, la ALB y en menor medida la A.G. Hay algunas variables que tienen la mayor parte de su representación en las dos o tres primeras componentes principales, como la ALB o la Spgt, mientras que otras como Age se encuentran en la quinta dimensión y otras repartidas en las primeras componentes y las últimas (A.G).

A parte del mapa de calor de la Figura 1.4, también podemos visualizar un gráfico de barras donde se nos da la calidad de representación de las variables en un número determinado de dimensiones. Por ejemplo, en la Figura 1.5, vemos cuatro gráficos de barras indicandonos la calidad de la representación total de las variables en las 2, 3, 5 o 6 primeras componentes principales.

themes <- theme(axis.title.y = element_text(size = 9))

cos1_2 <- fviz_cos2(pca, choice = "var", axes = 1:2, fill = "darkblue", color = "darkblue") + themes
cos1_3 <- fviz_cos2(pca, choice = "var", axes = 1:3, fill = "blue2", color = "blue2") + themes

cos1_5 <- fviz_cos2(pca, choice = "var", axes = 1:5, fill = "cadetblue", color = "cadetblue")  + themes
cos1_6 <- fviz_cos2(pca, choice = "var", axes = 1:6, fill = "lightblue", color = "lightblue") + themes

(cos1_2+cos1_3)/(cos1_5+cos1_6)
Figura 1.5. Calidad de represntacion de las variables en las 2, 3, 5 o 6 primeras componentes principales.

Figura 1.5. Calidad de represntacion de las variables en las 2, 3, 5 o 6 primeras componentes principales.

  • Algo que se ve claramente es que con las dos primeras dimensiones, hay una gran diferencia de calidad de representación entre la ALB (la mejor representada) y las otras variables, que mayoritariamanete se encuentran por encima de 0.5, salvo Age y Alkphos.

  • Si aumentamos el número de dimensiones y cogemos las tres primeras, lo que hemos considerado con el criterio del Kaiser (con un corte de 1), la calidad de la representación de la mayoría de las variables aumenta (salvo por A.G, Age y Alkphos).

  • Si aumentamos las dimensiones a las 5 primeras (que explican casi el 90% de la variablidad y hemos considerado con el criterio de Kaiser con un corte de 0.7), tenemos que casi todas las variables (salvo A.G) estan en un valor de cos2 mayor que 0.75 indicando una muy buena representación.

  • Finalmente, con las seis primeras componentes (que hemos escogido según el criterio gráfico ), todas las variables tienen una calidad de representación muy alta, con un cos2 cerca de 1.

Alternativamente, también se pueden visualizar los datos mediante las coordenadas por pares de dimensiones con la función fviz_pca_var(). Esta función no nos da exactamente la calidad de la representación, sinó que representa las variables usando las coordenadas para cada dimensión y visualizando las dimensiones de 2 en 2, de esta manera podemos observar las correlaciones entre variables en dichas dimensiones. Además, si consideramos un gradiente de colores en base al valor de cos2, podemos observar la calidad de la representación en las dimensiones observadas, tal y como se ha hecho en la Figura S1.2. Así mismo, este tipo de gráficos o los gráficos de barras del estilo de la Figura 1.5 pueden usarse para ver la contribución de cada variable a cada una de las dimensiones (funciones fviz_pca_contrib() y fviz_contrib()).

Cabe destacar que nos hemos centrado en la calidad de representación de las variables, pero de la misma forma, podríamos estudiar la calidad de la representación de los individuos si usaramos la funciones análogas a las usadas (i.e. fviz_cos2(pca, choice = "ind", ...), etc).




B)

Interpretar los dos primeros ejes principales mediante el gráfico de correlaciones con las variables originales.

El gráfico de correlaciones se obtiene con la función fviz_pca_var(), como se ha hecho en las figuras S1.2 y 1.6. Con este gráfico podemos ver las correlaciones entre las distintas variables en las dimensiones indicadas (en este caso PC1 y PC2):

  • Cuanto más cerca estan dos flechas en un cuadrante, más correlacionadas estan.
  • Las correlaciones positivas se encuentran en el mismo cuadrante.
  • Las correlaciones negativas se encuentran en cuadrantes opuestos.
  • La distancia entre el origen y la variable (punta de la flecha) indica la calidad de la variable en el mapa de factores (a más distancia, más calidad.).
fviz_pca_var(pca, col.var = "contrib", gradient.cols =c("#00AFBB", "#E7B800", "#FC4E07"))
Figura 1.6. Gráfico de correlaciones de las variables en las dos primeras componentes principales. El gradiente de color indica la contribución de cada variable a las dimensiones correspondientes.

Figura 1.6. Gráfico de correlaciones de las variables en las dos primeras componentes principales. El gradiente de color indica la contribución de cada variable a las dimensiones correspondientes.

En la Figura 1.6 vemos como hay algunas variables muy correlacionadas en las primeras dos dimensiones. Por ejemplo, TB y DB tienen una correlación casi perfecta en las componenentes principales 1 y 2, así como Sgot y Sgpt. En cambio, vemos coo la variable ALB está inversamente correalcionada con la variable Age, aunque vasándonos en su contribución y su calidad de representación, no podemos hacer predicciones sobre la variable Age en las componentes 1 y 2.




C)

Representar los pacientes en un gráfico de dos dimensiones con puntos distintos según el grupo y el género.

pca <- princomp(data_filt %>% dplyr::select(everything(), -Gender, -Selector), cor = T)

# mutate datafilt to create a new variable mixing gender and selector
group <- data_filt %>% mutate(SelGen = paste(Gender, "-", Selector))

fviz_pca_ind(pca,
             geom.ind = "point",
             fill.ind = group$SelGen,
             pointshape = 21, pointsize = 2,
             legend.title = list(fill = "Gender", color = "Selector"))
Figura 1.7. Gráfico en als dimensiones 1 (eje X) y 2 (eje Y) de los pacientes segun el grupo (Selector) y género (Gender).

Figura 1.7. Gráfico en als dimensiones 1 (eje X) y 2 (eje Y) de los pacientes segun el grupo (Selector) y género (Gender).




D)

¿Alguien se atreve con el gráfico en tres dimensiones?

El gráfico de el PCA en tres dimensiones se puede llevar a cabo mediante el paquete pca3d.

library(pca3d)

pc <- princomp(data_num, cor=TRUE, scores=TRUE)

gr <- data_filt$Gender
pca3d(pca = pc, group = gr, components = c(1,2,3), show.scale = T, show.plane = T, legend = "topleft", title = "3D-PCA plot on the first 3 PCs. Gender coloring.")
snapshotPCA3d(file="figs/PCA3d_gender.png")

gr <- data_filt$Selector
pca3d(pca = pc, group = gr, components = c(1,2,3), show.scale = T, show.plane = T, legend = "topleft", title = "3D-PCA plot on the first 3 PCs. Selector coloring.")
snapshotPCA3d(file="figs/PCA3d_selector.png")


gr <- data_filt %>% mutate(SelGen = paste(Gender, "-", Selector))
gr <- gr$SelGen
pca3d(pca = pc, group = gr, components = c(1,2,3), show.scale = T, show.plane = T, legend = "topleft", title = "3D-PCA plot on the first 3 PCs. Selector coloring.")
snapshotPCA3d(file="figs/PCA3d_SelGen.png")
knitr::include_graphics("figs/PCA3d_SelGen.png")
Figura 1.8. Gráfico 3D de los pacientes en las tres primeras componentes. Los puntos están coloresados según el género y el grupo de cada uno.

Figura 1.8. Gráfico 3D de los pacientes en las tres primeras componentes. Los puntos están coloresados según el género y el grupo de cada uno.

En la Figura S1.3 se puede observar un gráfico similar, pero usando la función scatterplot3d() del paquete con el mismo nombre.







Ejercicio 2

Una de las finalidades de disponer de marcadores para la enfermedad es determinar si ciertas variables se encuentran alteradas en los pacientes afectados. Para ello vamos a utilizar los datos del ejercicio 1, para investigar si hay diferencias entre los dos grupos.

A)

a. Realizar un boxplot múltiple que permita comparar visualmente todas las variables numéricas, separadas por controles y pacientes.

Para hacer un boxplot múltiple, primero de todo usaremos la dunción melt() del paquete reshape2 y después usaremos el paquete ggplot2 (funciones ggplot() y geom_boxplot()). Finalmente usaremos la función facet_wrap() para separar los valores de la variable Selector en dos paneles o, por otro lado, también podremos separar todas las variables en distintos paneles para hacer más fácil la comparación entre pacientes y controles por cada variable.

data_filt %>% 
  dplyr::mutate(Selector = Selector %>% recode("1" = "Pacients","2" = "Controls")) %>%
  
  melt() %>%
  
  ggplot(aes(x=Selector, y=value, fill=Selector)) + geom_boxplot(outlier.colour = "Gray50") +  
  
  facet_wrap(~variable, scales = "free_y") +
  
  stat_compare_means(method = "wilcox.test", paired = F,
                     label = "p.format", label.y.npc = 0.85, label.x.npc = 0.3) +

  theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 0.95),
        legend.position = "none") +
  
  ggtitle("Variables of our study separated in patients and control subjects.") +
  ylab("") + xlab("")
Figura 2.1. Gráficos de cajas de las variables numéricas presentes en los datos del ejercicio 1, separados según su condición ('Selector') de pacientes o controles. El p-valor mostrado ha sido obtenido mediante un test de wilcoxon (no pareado) e indica si las medias de pacientes y controles para la variable correspondiente son significativamente distintas o no.

Figura 2.1. Gráficos de cajas de las variables numéricas presentes en los datos del ejercicio 1, separados según su condición (‘Selector’) de pacientes o controles. El p-valor mostrado ha sido obtenido mediante un test de wilcoxon (no pareado) e indica si las medias de pacientes y controles para la variable correspondiente son significativamente distintas o no.

La Figura 2.1 nos permite comparar fácilmente muchas de las variables entre los pacientes (Selector = 1) y los controles (Selector = 2) de nuestro estudio:

  • La edad de los pacientes (mediana \(\simeq\) 46) es ligeramente mayor que los controles (mediana \(\simeq\) 41).

  • Aunque es dificil de ver a simple vista la DB y la TB de los pacientes (medianas: DB \(\simeq\) 0.5, TB \(\simeq\) 1.4) es ligeramente superior a la de los controles (medianas: DB \(\simeq\) 0.2, TB \(\simeq\) 0.8).

  • Por lo que respeta a las variables TP, ALB y A.G, se ve una ligera diferencia, siendo sus valores mayores en el caso del grupo control (medianas: TP \(\simeq\) 6.6, ALB \(\simeq\) 3.4, A.G \(\simeq\) 1), que en el grupo de los pacientes (medianas: TP \(\simeq\) 6.55, ALB \(\simeq\) 3, A.G \(\simeq\) 0.9). Sin embargo, esta diferencia no es significativa para la variable TP.

  • Finalmente están las variables Alkphos, Sgpt y Sgot, las medias són significativamente diferentes y con unos valores mayores en el caso de los pacientes (medianas: Alkphos \(\simeq\) 229, Sgpt \(\simeq\) 51, Sgot \(\simeq\) 53) que de los sujetos del grupo control (medianas: Alkphos \(\simeq\) 187, Sgpt \(\simeq\) 28, Sgot \(\simeq\) 29).

En la Figura S2.1 se puede ver un gráfico con los mismos datos separados en dos paneles, uno para pacientes y otro para controles.







B)

Si suponemos que se verifican todas las condiciones, realizar un MANOVA de un factor, utilizando la variable Selector, para decidir si hay diferencias en las variables marcadores

Antes de realizar el test MANOVA debemos hacer varias suposciciones:

  • Normalidad: esta se puede comprobar mediante el test de Shapiro-Wilk multivariante, con la función mshapiro_test().
mshapiro_test(data = data_num)
## # A tibble: 1 x 2
##   statistic  p.value
##       <dbl>    <dbl>
## 1     0.124 2.55e-45

  • Homogeneidad de las varianzas: esto se puede llevar a cabo con un test de Levene (varianzas) o con un test de Box (covarianzas). Más adelante veremos que las varianzas y covarianzas entre grupos no son homogeneas.

  • …

Como hemos visto, los datos no cumplen una distribución normal, ya que el p-valor obtenido del test de Shapiro-Wilk, es mucho menor que el nivel de significancia de 0.05, con lo que no podemos asumir la hipòtesis nula de la distribución normal. Así pues, tendríamos que usar un test no paramétrico homólogo al MANOVA.

Sin embargo, asumiendo que todas las condiciones mencionadas se cumplen, para hacer el test MANOVA en un factor (Selector) podemos usar la función manova() que admite como argumentos los mismos que la función aov para el ANOVA.

attach(data_filt)
manova <- manova(cbind(A.G, ALB, TP, Sgot, Sgpt, Alkphos, DB, TB, Age) ~ Selector, data = data_filt)
detach()
summary(manova)
##            Df  Pillai approx F num Df den Df    Pr(>F)    
## Selector    1 0.11654     8.34      9    569 9.984e-12 ***
## Residuals 577                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Segun el test MANOVA que hemos realizado, sí que hay diferencias significativas en las variables marcadores.

summary.aov(manova) 
##  Response A.G :
##              Df Sum Sq Mean Sq F value    Pr(>F)    
## Selector      1  1.571 1.57107  15.775 8.037e-05 ***
## Residuals   577 57.465 0.09959                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response ALB :
##              Df Sum Sq Mean Sq F value    Pr(>F)    
## Selector      1   9.31  9.3118  15.114 0.0001129 ***
## Residuals   577 355.48  0.6161                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response TP :
##              Df Sum Sq Mean Sq F value Pr(>F)
## Selector      1   0.77 0.76831  0.6527 0.4195
## Residuals   577 679.22 1.17715               
## 
##  Response Sgot :
##              Df   Sum Sq Mean Sq F value    Pr(>F)    
## Selector      1  1119477 1119477  13.616 0.0002455 ***
## Residuals   577 47440061   82218                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response Sgpt :
##              Df   Sum Sq Mean Sq F value    Pr(>F)    
## Selector      1   516055  516055  15.772 8.049e-05 ***
## Residuals   577 18879287   32720                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response Alkphos :
##              Df   Sum Sq Mean Sq F value    Pr(>F)    
## Selector      1  1152846 1152846  20.075 8.982e-06 ***
## Residuals   577 33135491   57427                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response DB :
##              Df Sum Sq Mean Sq F value    Pr(>F)    
## Selector      1  278.1 278.087  37.255 1.903e-09 ***
## Residuals   577 4307.0   7.464                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response TB :
##              Df  Sum Sq Mean Sq F value    Pr(>F)    
## Selector      1  1087.2 1087.15  29.408 8.633e-08 ***
## Residuals   577 21330.3   36.97                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response Age :
##              Df Sum Sq Mean Sq F value  Pr(>F)   
## Selector      1   2697 2697.09  10.416 0.00132 **
## Residuals   577 149401  258.93                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Concretamente, con un nivel de significancia de 0.05, todas las variables numéricas salvo TP muestran diferencias significativas entre los dos grupos.







C)

Como además tenemos el factor Gender podemos considerar un modelo general con todos los factores (sin interacciones) y el modelo más simple con solo el factor principal y contrastar con la función anova() si el modelo simple explica la variación de los datos.

Para hacer el modelo general con los dos factores (Selector y Gender) usaremos la función manova(), como hicimos en el apartado B, donde generamos el modelo simple con un solo factor (Selector). Luego para comparar los dos modelos, podemos usar la función summary() o la función anova() y ver que factores explican significativamente la variación de los datos.

attach(data_filt)
manova_general <- manova(cbind(A.G, ALB, TP, Sgot, Sgpt, Alkphos, DB, TB, Age) ~ Selector+Gender, data = data_filt)
detach()

Con el model simple de un factor:

anova(manova)
## Analysis of Variance Table
## 
##              Df  Pillai approx F num Df den Df    Pr(>F)    
## (Intercept)   1 0.98774   5095.5      9    569 < 2.2e-16 ***
## Selector      1 0.11654      8.3      9    569 9.984e-12 ***
## Residuals   577                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Vemos como la variable Selector explica significativamente la variación de los datos.

Con el model general con los dos factores:

anova(manova_general)
## Analysis of Variance Table
## 
##              Df  Pillai approx F num Df den Df    Pr(>F)    
## (Intercept)   1 0.98775   5086.9      9    568 < 2.2e-16 ***
## Selector      1 0.11717      8.4      9    568 8.788e-12 ***
## Gender        1 0.02898      1.9      9    568   0.05186 .  
## Residuals   576                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Vemos como la variable Selector explica significativamente la variación de los datos, mientras que la variable Gender no es significativa (usando un nivel de significancia de 0.05, pero sí si usamos de 0.1).

Así pues, como en el modelo general el único factor que explica la variación de los datos es Selector, podemos concluir que con el modelo simple tenemos suficiente para explicar dicha variación.







D)

Comparar las matrices de covarianzas de los dos grupos del factor Selector con el test de la razón de verosimilitudes y con el test M de Box

Las matrices de covarianzas se pueden encontrar en la Figura S2.1 para Selector = 1 y en la Figura S2.2 para Selector = 2.

  • Test de razón de verosimilitudes:

  • Test M de Box: el test M de Box, o prueba de box permite comparar las matrices de covarianzas entre dos grupos para ver si estas son iguales. Es un test muy sensible y debido a esto, se recomienda un límite de significancia bastante bajo (i.e. 0.001). También hay que tener en cuenta, que en caso de falta de normalidad de los datos, el resultado del test de Box puede ser significativo debido a esta falta de normalidad y no a la diferencia entre las matrices de covarianzas. Para llevar a cabo el test, se puede usar la función boxM() del paquete heplots.

boxM_res <- heplots::boxM(data_num, data_filt$Selector)
summary(boxM_res)
## Summary for Box's M-test of Equality of Covariance Matrices
## 
## Chi-Sq:   2496.705 
## df:   45 
## p-value: < 2.2e-16 
## 
## log of Covariance determinants:
##        1        2   pooled 
## 38.00520 16.17060 36.22022 
## 
## Eigenvalues:
##              1            2       pooled
## 1 1.488917e+05 2.023218e+04 1.074514e+05
## 2 6.973537e+04 1.505127e+03 5.519672e+04
## 3 1.346050e+04 2.889142e+02 9.726111e+03
## 4 2.378558e+02 2.647388e+02 2.532873e+02
## 5 5.525106e+01 1.527346e+00 3.979737e+01
## 6 2.193000e+00 7.650083e-01 1.573901e+00
## 7 1.542026e+00 1.544093e-01 1.548583e+00
## 8 2.083380e-01 5.788200e-03 1.973065e-01
## 9 2.474629e-02 4.332887e-03 1.921571e-02
## 
## Statistics based on eigenvalues:
##                      1            2       pooled
## product   3.202195e+16 1.053903e+07 5.373311e+15
## sum       2.323846e+05 2.229342e+04 1.726707e+05
## precision 2.158108e-02 2.427158e-03 1.711762e-02
## max       1.488917e+05 2.023218e+04 1.074514e+05

Como se puede observar, con la prueba de Box, obtenemos que las matrices de covarianzas de los dos grupos son significativamente distintas. Aunque bien podría deberse a la falta de normalidad de los datos.







E)

Otra opción para el apartado anterior habría sido comparar la variabilidad con un test de Levene múltiple. ¿Este test contrasta lo mismo que los anteriores?

El test de Levene es un test diseñado para el análisis de varianzas. Está diseñado para análisis univariante de la varianza (ANOVA), aunque se puede extrapolar a multivariante, pero con este contrastamos la homogeneidad de cada variable (igualdad de varianzas) entre los grupos, mientras con los dos tests efectuados en el apartado D, miramos la varianza y las covarianzas de todas las variables entre los grupos estudiados (Selector = 1, Selector = 2).

library(car)
attach(data_filt)
levene_test(data = data_filt, formula = A.G+ALB+TP+Sgot+Sgpt+Alkphos+DB+TB+Age ~ Selector)
## # A tibble: 1 x 4
##     df1   df2 statistic          p
##   <int> <int>     <dbl>      <dbl>
## 1     1   577      23.5 0.00000158
detach()







F)

Los experimentadores dudan de la normalidad de las variables observadas. Estudiar el caso.

f1. Estudiar gráficamente y con algún test la normalidad univariante y multivariante (si es posible) en cada uno de los grupos

Gràficamente, podemos usar la función mvn() para hacer los gráficos univariantes (i.e. histogramas) y multivariante (i.e. qqplot):

mvn_res <- MVN::mvn(data = data_num_selector1, univariatePlot = "histogram")
Figura 2.2a. Histogramas de para cada variable en los datos de pacientes (Selector = 1).

Figura 2.2a. Histogramas de para cada variable en los datos de pacientes (Selector = 1). 

mvn_res <- MVN::mvn(data = data_num_selector2, univariatePlot = "histogram")
Figura 2.2b. Histogramas de para cada variable en los datos de sujetos control (Selector = 2).

Figura 2.2b. Histogramas de para cada variable en los datos de sujetos control (Selector = 2). 

mvn_res <- MVN::mvn(data = data_num_selector1, multivariatePlot = "qq")
Figura 2.3a. QQplot de las distancias de Mahalanobis en los datos de pacientes (Selector = 1)

Figura 2.3a. QQplot de las distancias de Mahalanobis en los datos de pacientes (Selector = 1) 

mvn_res <- MVN::mvn(data = data_num_selector2, multivariatePlot = "qq")
Figura 2.3b. QQplot de las distancias de Mahalanobis en los datos de sujetos control (Selector = 2)

Figura 2.3b. QQplot de las distancias de Mahalanobis en los datos de sujetos control (Selector = 2)

A nivel gráfico, observando las Figuras 2.2a y 2.2b se podría intuir una normalidad univariante para las variables Age y TP en ambos grupos, mientras que ALB parece tener una distribución similar a la normal en el sólo grupo de los pacientes y las otras variable no se parecen para nada a una distribución normal, pues són altamente asimétricas en ambos grupos.

En cuanto a los gráficos (qqplots) multivariantes (Figuras 2.3a y 2.3b), estos nos dicen que los datos no siguen una normalidad multivariante.

Sin embargo, no podemos fiarnos solo de los métodos gráficos y se tienen que usar otros test estadísticos, como Shapiro-Wilk, que nos servirá tanto a nivel univariante como multivariante. * Shapiro-Wilk univariante:

shapiro_selector1 <- shapiro_test(data_num_selector1, vars = colnames(data_num_selector1))
shapiro_selector2 <- shapiro_test(data_num_selector2, vars = colnames(data_num_selector1))

paste("Resultado del test Shapiro-Wilk para el grupo Selector 1:")
## [1] "Resultado del test Shapiro-Wilk para el grupo Selector 1:"
shapiro_selector1
## # A tibble: 9 x 3
##   variable statistic        p
##   <chr>        <dbl>    <dbl>
## 1 A.G          0.928 3.31e-13
## 2 Age          0.990 6.28e- 3
## 3 ALB          0.993 4.59e- 2
## 4 Alkphos      0.624 4.29e-29
## 5 DB           0.607 1.39e-29
## 6 Sgot         0.320 2.16e-36
## 7 Sgpt         0.374 2.61e-35
## 8 TB           0.529 9.54e-32
## 9 TP           0.991 1.32e- 2
paste("Como se puede ver, todas las variables muestran un p-valor inferior a 0.05, por lo cual no podemos considerar que dichas variables sigan una distribución normal en el caso del grupo Selector 1.")
## [1] "Como se puede ver, todas las variables muestran un p-valor inferior a 0.05, por lo cual no podemos considerar que dichas variables sigan una distribución normal en el caso del grupo Selector 1."
paste("")
## [1] ""
paste("Resultado del test Shapiro-Wilk para el grupo Selector 2:")
## [1] "Resultado del test Shapiro-Wilk para el grupo Selector 2:"
shapiro_selector2
## # A tibble: 9 x 3
##   variable statistic        p
##   <chr>        <dbl>    <dbl>
## 1 A.G          0.969 9.15e- 4
## 2 Age          0.983 4.33e- 2
## 3 ALB          0.980 1.72e- 2
## 4 Alkphos      0.495 1.42e-21
## 5 DB           0.506 2.33e-21
## 6 Sgot         0.635 1.28e-18
## 7 Sgpt         0.686 2.54e-17
## 8 TB           0.500 1.75e-21
## 9 TP           0.989 2.21e- 1
paste("Como se puede ver, todas las variables muestran menos TP un p-valor inferior a 0.05, por lo cual no podemos considerar que dichas variables sigan una distribución normal en el caso del grupo Selector 2.")
## [1] "Como se puede ver, todas las variables muestran menos TP un p-valor inferior a 0.05, por lo cual no podemos considerar que dichas variables sigan una distribución normal en el caso del grupo Selector 2."
  • Shapiro-Wilk multivariante:
mshapiro_selector1 <- mshapiro_test(data = data_num_selector1)
mshapiro_selector2 <- mshapiro_test(data = data_num_selector2)

paste("El p-valor del test multivariante de Shapiro-Wilk en el grupo Selector 1 es", mshapiro_selector1$p.value, "por lo que no podemos considerar que este grupo siga una distribución normal multivariante.")
## [1] "El p-valor del test multivariante de Shapiro-Wilk en el grupo Selector 1 es 2.00419808624749e-39 por lo que no podemos considerar que este grupo siga una distribución normal multivariante."
paste("El p-valor del test multivariante de Shapiro-Wilk en el grupo Selector 2 es", mshapiro_selector2$p.value, "por lo que no podemos considerar que este grupo siga una distribución normal multivariante.")
## [1] "El p-valor del test multivariante de Shapiro-Wilk en el grupo Selector 2 es 1.16805887608274e-21 por lo que no podemos considerar que este grupo siga una distribución normal multivariante."

f2. La normalidad multivariante se puede estudiar con la asimetría y la kurtosis multivariante de Mardia gracias a la función mvn() del paquete MVN. La misma función permite el estudio univariante y también hace los gráficos.

mvn_res <- MVN::mvn(data = data_num)
mvn_res$multivariateNormality %>% knitr::kable(caption = "Tabla 2.1. Resultado de la función 'mvn()' para la normalidad multivariante.")
Tabla 2.1. Resultado de la función ‘mvn()’ para la normalidad multivariante.
Test Statistic p value Result
Mardia Skewness 71910.4739009141 0 NO
Mardia Kurtosis 788.563271538733 0 NO
MVN NA NA NO 
mvn_res$univariateNormality %>% knitr::kable(caption = "Tabla 2.2. Resultado de la función 'mvn()' para la normalidad univariante.")
Tabla 2.2. Resultado de la función ‘mvn()’ para la normalidad univariante.
Test Variable Statistic p value Normality
Shapiro-Wilk Age 0.9920 0.0033 NO
Shapiro-Wilk TB 0.4613 <0.001 NO
Shapiro-Wilk DB 0.5312 <0.001 NO
Shapiro-Wilk Alkphos 0.5859 <0.001 NO
Shapiro-Wilk Sgpt 0.3286 <0.001 NO
Shapiro-Wilk Sgot 0.2811 <0.001 NO
Shapiro-Wilk TP 0.9918 0.0029 NO
Shapiro-Wilk ALB 0.9925 0.0053 NO
Shapiro-Wilk A.G 0.9465 <0.001 NO

f3. Acompañar el test con un gráfico qq-plot de ajuste a la distribución ji-cuadrado de las distancias de Mahalanobis (clásicas y robustas) de los datos a la media en cada grupo.

d2 <- Moutlier(data_num, quantile = 0.975, plot = F)
n <- nrow(data_num)
p <- ncol(data_num)
quantiles <- qchisq((1:n)/(n+1),p)
plot(quantiles, sort(d2$md), ylab = "Ordered classical Mahalanobis distances",xlab = "Chi-square quantile")
Figura 2.4. QQplot de ajuste a la distribución ji-cuadrado de las distancias de Mahalanobis clàsicas.

Figura 2.4. QQplot de ajuste a la distribución ji-cuadrado de las distancias de Mahalanobis clàsicas. 

plot(quantiles, sort(d2$rd), ylab = "Ordered robust Mahalanobis distances", xlab = "Chi-square quantile")
Figura 2.5. QQplot de ajuste a la distribución ji-cuadrado de las distancias de Mahalanobis robustas.

Figura 2.5. QQplot de ajuste a la distribución ji-cuadrado de las distancias de Mahalanobis robustas.

f4. Comentar el resultado en cuanto a los tests y en cuanto a la presencia de datos atípicos.

f5. En caso de duda sobre la normalidad multivariante de los datos, proponer y calcular un método multivariante para contrastar si hay diferencias de variabilidad entre los grupos.







G)

Según cual haya sido el resultado de los apartados anteriores y dado que son únicamente dos grupos, tal vez se deba utilizar un test de permutaciones con la T2 de Hotelling (Estudiar la función hotelling.test() del paquete Hotelling)







Material suplementario

Tablas

cor_list_1d %>% knitr::kable(align = "c", caption = "Tabla S1.1. Lista de correlaciones dos a dos, ordenadas de mayor a menor correlación (en valor absoluto). Los valores de 1, indicando que una variable se enfrenta a ella misma, son mantenidos.")
Tabla S1.1. Lista de correlaciones dos a dos, ordenadas de mayor a menor correlación (en valor absoluto). Los valores de 1, indicando que una variable se enfrenta a ella misma, son mantenidos.
Variable 1 Variable 2 Pearson’s correlation Absolute Pearson’s
Age Age 1.0000000 1.0000000
TB TB 1.0000000 1.0000000
DB DB 1.0000000 1.0000000
Alkphos Alkphos 1.0000000 1.0000000
Sgpt Sgpt 1.0000000 1.0000000
Sgot Sgot 1.0000000 1.0000000
TP TP 1.0000000 1.0000000
ALB ALB 1.0000000 1.0000000
A.G A.G 1.0000000 1.0000000
DB TB 0.8744810 0.8744810
TB DB 0.8744810 0.8744810
Sgot Sgpt 0.7918621 0.7918621
Sgpt Sgot 0.7918621 0.7918621
ALB TP 0.7831122 0.7831122
TP ALB 0.7831122 0.7831122
A.G ALB 0.6896323 0.6896323
ALB A.G 0.6896323 0.6896323
ALB Age -0.2642109 0.2642109
Age ALB -0.2642109 0.2642109
Sgot DB 0.2570224 0.2570224
DB Sgot 0.2570224 0.2570224
Sgot TB 0.2373231 0.2373231
TB Sgot 0.2373231 0.2373231
A.G TP 0.2348872 0.2348872
TP A.G 0.2348872 0.2348872
A.G Alkphos -0.2341665 0.2341665
Alkphos A.G -0.2341665 0.2341665
Alkphos DB 0.2340076 0.2340076
DB Alkphos 0.2340076 0.2340076
Sgpt DB 0.2331801 0.2331801
DB Sgpt 0.2331801 0.2331801
ALB DB -0.2284092 0.2284092
DB ALB -0.2284092 0.2284092
ALB TB -0.2220866 0.2220866
TB ALB -0.2220866 0.2220866
A.G Age -0.2164083 0.2164083
Age A.G -0.2164083 0.2164083
Sgpt TB 0.2133755 0.2133755
TB Sgpt 0.2133755 0.2133755
A.G TB -0.2062672 0.2062672
TB A.G -0.2062672 0.2062672
Alkphos TB 0.2057392 0.2057392
TB Alkphos 0.2057392 0.2057392
A.G DB -0.2001247 0.2001247
DB A.G -0.2001247 0.2001247
TP Age -0.1862481 0.1862481
Age TP -0.1862481 0.1862481
Sgot Alkphos 0.1665800 0.1665800
Alkphos Sgot 0.1665800 0.1665800
ALB Alkphos -0.1634186 0.1634186
Alkphos ALB -0.1634186 0.1634186
Sgpt Alkphos 0.1247767 0.1247767
Alkphos Sgpt 0.1247767 0.1247767
Sgpt Age -0.0877992 0.0877992
Age Sgpt -0.0877992 0.0877992
ALB Sgot -0.0849146 0.0849146
Sgot ALB -0.0849146 0.0849146
Alkphos Age 0.0788784 0.0788784
Age Alkphos 0.0788784 0.0788784
A.G Sgot -0.0700398 0.0700398
Sgot A.G -0.0700398 0.0700398
TP Sgpt -0.0424321 0.0424321
Sgpt TP -0.0424321 0.0424321
ALB Sgpt -0.0286575 0.0286575
Sgpt ALB -0.0286575 0.0286575
TP Alkphos -0.0270620 0.0270620
Alkphos TP -0.0270620 0.0270620
TP Sgot -0.0257510 0.0257510
Sgot TP -0.0257510 0.0257510
Sgot Age -0.0204989 0.0204989
Age Sgot -0.0204989 0.0204989
TB Age 0.0110004 0.0110004
Age TB 0.0110004 0.0110004
TP TB -0.0079059 0.0079059
TB TP -0.0079059 0.0079059
DB Age 0.0067843 0.0067843
Age DB 0.0067843 0.0067843
A.G Sgpt -0.0023750 0.0023750
Sgpt A.G -0.0023750 0.0023750
TP DB 0.0000327 0.0000327
DB TP 0.0000327 0.0000327 
par_cor_1e_list <- par_cor_1e$estimate %>% 
  reshape2::melt() %>%
  dplyr::mutate(abs_pearson = sqrt(value^2)) %>%
  dplyr::arrange(desc(abs_pearson)) %>%
  magrittr::set_colnames(c("Variable 1", "Variable 2", "Pearson's correlation", "Absolute Pearson's"))

par_cor_1e_list %>% knitr::kable(align = "c", caption = "Tabla 1.2. Lista de correlaciones parciales dos a dos, ordenadas de mayor a menor correlación (en valor absoluto). Los valores de 1, indicando que una variable se enfrenta a ella misma, son mantenidos.")
Tabla 1.2. Lista de correlaciones parciales dos a dos, ordenadas de mayor a menor correlación (en valor absoluto). Los valores de 1, indicando que una variable se enfrenta a ella misma, son mantenidos.
Variable 1 Variable 2 Pearson’s correlation Absolute Pearson’s
Age Age 1.0000000 1.0000000
TB TB 1.0000000 1.0000000
DB DB 1.0000000 1.0000000
Alkphos Alkphos 1.0000000 1.0000000
Sgpt Sgpt 1.0000000 1.0000000
Sgot Sgot 1.0000000 1.0000000
TP TP 1.0000000 1.0000000
ALB ALB 1.0000000 1.0000000
A.G A.G 1.0000000 1.0000000
ALB TP 0.8971925 0.8971925
TP ALB 0.8971925 0.8971925
DB TB 0.8401910 0.8401910
TB DB 0.8401910 0.8401910
A.G ALB 0.8246399 0.8246399
ALB A.G 0.8246399 0.8246399
Sgot Sgpt 0.7829975 0.7829975
Sgpt Sgot 0.7829975 0.7829975
A.G TP -0.6865110 0.6865110
TP A.G -0.6865110 0.6865110
DB ALB -0.1910744 0.1910744
ALB DB -0.1910744 0.1910744
TP DB 0.1823035 0.1823035
DB TP 0.1823035 0.1823035
TP Sgpt -0.1656768 0.1656768
Sgpt TP -0.1656768 0.1656768
ALB Sgpt 0.1469404 0.1469404
Sgpt ALB 0.1469404 0.1469404
A.G DB 0.1412185 0.1412185
DB A.G 0.1412185 0.1412185
TP Sgot 0.1234029 0.1234029
Sgot TP 0.1234029 0.1234029
ALB Sgot -0.1171630 0.1171630
Sgot ALB -0.1171630 0.1171630
Alkphos A.G -0.1165394 0.1165394
A.G Alkphos -0.1165394 0.1165394
Alkphos DB 0.0971593 0.0971593
DB Alkphos 0.0971593 0.0971593
Age Sgpt -0.0957684 0.0957684
Sgpt Age -0.0957684 0.0957684
ALB Age -0.0820090 0.0820090
Age ALB -0.0820090 0.0820090
A.G Sgpt -0.0715968 0.0715968
Sgpt A.G -0.0715968 0.0715968
Alkphos Sgot 0.0700443 0.0700443
Sgot Alkphos 0.0700443 0.0700443
A.G Sgot 0.0591914 0.0591914
Sgot A.G 0.0591914 0.0591914
Sgpt DB 0.0579524 0.0579524
DB Sgpt 0.0579524 0.0579524
Sgot Age 0.0523929 0.0523929
Age Sgot 0.0523929 0.0523929
Age Alkphos 0.0461958 0.0461958
Alkphos Age 0.0461958 0.0461958
A.G TB -0.0293029 0.0293029
TB A.G -0.0293029 0.0293029
Age DB -0.0248805 0.0248805
DB Age -0.0248805 0.0248805
A.G Age -0.0228778 0.0228778
Age A.G -0.0228778 0.0228778
Alkphos ALB 0.0146719 0.0146719
ALB Alkphos 0.0146719 0.0146719
Alkphos TB -0.0135221 0.0135221
TB Alkphos -0.0135221 0.0135221
Sgot TB 0.0118750 0.0118750
TB Sgot 0.0118750 0.0118750
DB Sgot 0.0105666 0.0105666
Sgot DB 0.0105666 0.0105666
ALB TB -0.0087443 0.0087443
TB ALB -0.0087443 0.0087443
TB TP 0.0079615 0.0079615
TP TB 0.0079615 0.0079615
TP Age 0.0075072 0.0075072
Age TP 0.0075072 0.0075072
TB Sgpt 0.0061201 0.0061201
Sgpt TB 0.0061201 0.0061201
TB Age -0.0029639 0.0029639
Age TB -0.0029639 0.0029639
Sgpt Alkphos 0.0010541 0.0010541
Alkphos Sgpt 0.0010541 0.0010541
Alkphos TP -0.0002238 0.0002238
TP Alkphos -0.0002238 0.0002238 
cov_mat_selector1 %>% knitr::kable(align = "c", caption = "Tabla S2.1. Matriz de covarianzas del grupo de pacientes ('Selector' = 1).")
Tabla S2.1. Matriz de covarianzas del grupo de pacientes (‘Selector’ = 1).
Age TB DB Alkphos Sgpt Sgot TP ALB A.G
Age 246.187178 -2.6722181 -1.7319577 376.866407 -436.8823385 -243.458680 -3.1780187 -3.4973541 -1.3308976
TB -2.672218 51.2423524 19.9290225 318.895747 279.9389462 509.423992 0.0609333 -1.2364565 -0.4557098
DB -1.731958 19.9290225 10.3203765 164.754420 137.6199711 247.996179 0.0670728 -0.5629437 -0.1916362
Alkphos 376.866407 318.8957469 164.7544198 72277.363091 5390.7891006 13157.741868 -5.9794926 -29.9190722 -18.3494859
Sgpt -436.882339 279.9389462 137.6199711 5390.789101 45461.4641834 56759.007053 -10.6640541 -0.8768935 1.9592638
Sgot -243.458680 509.4239920 247.9961785 13157.741868 56759.0070534 114336.114960 -7.4678598 -18.6288346 -5.9487430
TP -3.178019 0.0609333 0.0670728 -5.979493 -10.6640541 -7.467860 1.2040283 0.6576987 0.0720592
ALB -3.497354 -1.2364565 -0.5629437 -29.919072 -0.8768935 -18.628835 0.6576987 0.6200131 0.1702830
A.G -1.330898 -0.4557098 -0.1916362 -18.349486 1.9592638 -5.948743 0.0720592 0.1702830 0.1063755 
cov_mat_selector2  %>% knitr::kable(align = "c", caption = "Tabla S2.2. Matriz de covarianzas del grupo control ('Selector' = 2).")
Tabla S2.2. Matriz de covarianzas del grupo control (‘Selector’ = 2).
Age TB DB Alkphos Sgpt Sgot TP ALB A.G
Age 291.0133038 0.2049335 0.1732816 -190.695676 -46.7998891 -61.6452328 -3.2686807 -2.2266075 -0.2056375
TB 0.2049335 1.0193178 0.5197982 80.923976 8.5153104 14.2240798 -0.1654361 -0.1452531 -0.0472918
DB 0.1732816 0.5197982 0.2724257 41.683603 4.3908647 7.3869401 -0.0794290 -0.0732705 -0.0248308
Alkphos -190.6956763 80.9239763 41.6836031 20030.119586 1341.8261641 1384.5409091 -4.3997044 -16.1198263 -9.8257443
Sgpt -46.7998891 8.5153104 4.3908647 1341.826164 632.3328160 610.3452328 0.9814856 0.8766075 0.0663326
Sgot -61.6452328 14.2240798 7.3869401 1384.540909 610.3452328 1336.8645233 -4.0711197 -2.3125831 0.2007528
TP -3.2686807 -0.1654361 -0.0794290 -4.399704 0.9814856 -4.0711197 1.1094752 0.7056338 0.0987973
ALB -2.2266075 -0.1452531 -0.0732705 -16.119826 0.8766075 -2.3125831 0.7056338 0.6061826 0.1649558
A.G -0.2056375 -0.0472918 -0.0248308 -9.825744 0.0663326 0.2007528 0.0987973 0.1649558 0.0825138

Figuras

corrplot::corrplot(par_cor_1e$estimate, type = "lower", method = "color", 
                   addgrid.col = "gray50", addCoef.col = "black", 
                   title = "Partial correlation between all numeric variables one-vs-one.")
Figura S1.1. Mapa de calor con las correlaciones parciales entre todas las variables por pares. Las correlaciones parciales se han obtenido con la función 'pcor()' del paquete 'ppcor' y el gráfico se ha hecho con la función 'corrplot' del paquete 'corrplot'.

Figura S1.1. Mapa de calor con las correlaciones parciales entre todas las variables por pares. Las correlaciones parciales se han obtenido con la función ‘pcor()’ del paquete ‘ppcor’ y el gráfico se ha hecho con la función ‘corrplot’ del paquete ‘corrplot’. 

pca <- princomp(x = data_num, cor = T)

fviz_pca_var(pca, col.var = "cos2", 
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel = TRUE, axes = c(1,2), title="PC1 - PC2", labelsize = 4) +

fviz_pca_var(pca, col.var = "cos2", 
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel = TRUE, axes = c(1,3), title="PC1 - PC3", labelsize = 4) +

fviz_pca_var(pca, col.var = "cos2", 
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel = TRUE, axes = c(2,3), title="PC2 - PC3", labelsize = 4)
Figura S1.2.

Figura S1.2. 

library(scatterplot3d)
scatterplot3d(pc$scores[,1], pc$scores[,2], pc$scores[,3],
              xlab="PC1", ylab="PC2", zlab="PC3",
              color = ifelse(data_filt$Selector == "1", "green", "yellow"),
              pch = ifelse(data_filt$Gender == "Male", 1, 15), 
              main = "3D PCA with the first 3 components colored by selector and gender.")
Figura S1.3. PCA en 3D con las primeras tres componentes principales, con los individuos coloreados por la variable selector y con distintas formas para la variable Gender.

Figura S1.3. PCA en 3D con las primeras tres componentes principales, con los individuos coloreados por la variable selector y con distintas formas para la variable Gender. 

data_filt %>% 
  dplyr::mutate(Selector = Selector %>% dplyr::recode("1" = "Patients",
                                                      "2" = "Controls")) %>%
  
  melt() %>%
  
  ggplot(aes(x=variable, y=value, fill=variable)) + geom_boxplot(outlier.colour = "Gray50") +
  
  coord_cartesian(ylim = c(0,150)) +
  
  facet_wrap(~Selector) +
  
  theme(axis.text.x = element_text(angle = 30, hjust = 1, vjust = 0.95),
        legend.position = "none") +
  
  ggtitle("Variables of our study separated in patients and control subjects.") +
  ylab("") + xlab("")
Figura S2.1. Gráficos de cajas de las variables numéricas presentes en los datos del ejercicio 1, separados en dos paneles según su condición ('Selector') de pacienteso controles. El eje Y ha sido limitado para una mejor visualización de las cajas a cambio de no visualizar algunos outliers de algunas variables.

Figura S2.1. Gráficos de cajas de las variables numéricas presentes en los datos del ejercicio 1, separados en dos paneles según su condición (‘Selector’) de pacienteso controles. El eje Y ha sido limitado para una mejor visualización de las cajas a cambio de no visualizar algunos outliers de algunas variables.

Bibliografía consultada

  • Everitt and Hothorn (2011). An introduction to applied multivariate analysis with R. Springer.

  • C.M. Cuadras (2018). Nuevos métodos de análisis multivariante. CMC Editors.

  • Elsieo Martínez H. Tratamiento matricial de los datos multivariantes. Disponible online.

  • Kassambara (2017). Practical guide to principal component analysis in R. STHDA. Disponible online.

  • Pryiank Goyal (2020). PCA. Rpubs. Disponible online.